# 1. Define Population Distribution (e.g., skewed population)
set.seed(123)
population <- rgamma(100000, shape = 2, scale = 2)

# 2. Take a Sample Distribution (e.g., 100 individuals)
sample_data <- sample(population, 100)

set.seed(1)
OBV <- 1:10
Dist1 <- NULL
Dist9 <- NULL
Dist16 <- NULL
Dist25 <- NULL
Dist36 <- NULL
count = 100
while(count > 0){Dist1 <- c(Dist1,sample(OBV, 1, replace = TRUE)); count <- count - 1}
count = 100
while(count > 0){Dist9 <- c(Dist9,mean(sample(OBV, 9,replace = TRUE) ) ); count <- count - 1}
count = 100
while(count > 0){Dist16 <- c(Dist16,mean(sample(OBV, 16,replace = TRUE) ) ); count <- count - 1}
count = 100
while(count > 0){Dist25 <- c(Dist25,mean(sample(OBV, 25,replace = TRUE) ) ); count <- count - 1}
count = 100
while(count > 0){Dist36 <- c(Dist36,mean(sample(OBV, 36,replace = TRUE) ) ); count <- count - 1}
Dist.df <- data.frame(Size = factor(rep(c(1,9,16,25,36), each=100)), Sample_Means = c(Dist1, Dist9, Dist16, Dist25, Dist36) )
ggplot(Dist.df, aes(Sample_Means, fill = Size)) + geom_histogram() + facet_grid(. ~ Size)

set.seed(42)  # for reproducibility

# Number of coin flips
n_flips <- 10000

# Simulate rolling a fair 6-sided die
flips <- sample(c("H", "T"), size = n_flips, replace = TRUE)

# Compute cumulative mean of rolling a '1'
cumulative_mean <- cumsum(flips == "H") / (1:n_flips)

# Plot convergence
plot(1:n_flips, cumulative_mean, type = "l", col = "blue", lwd = 2,
     xlab = "Number of Flips", ylab = "Proportion of Heads",
     main = "Law of Large Numbers: Convergence to 1/2")
abline(h = 1/2, col = "red", lty = 2, lwd = 2)  # Reference line at 1/2

set.seed(42)  # For reproducibility

# Number of dice rolls
n_flips <- 10000

# Simulate rolling a fair 6-sided die
flips <- sample(c("H", "T"), size = n_flips, replace = TRUE)

# Prior: Beta(1,5) (weak prior belief about p = 1/6)
alpha <- 1  # prior successes (rolling a 1)
beta <- 2   # prior failures (rolling 2-6)

# Store posterior mean estimates
posterior_means <- numeric(n_flips)

# Bayesian updating
for (i in 1:n_flips) {
  alpha <- alpha + (flips[i] == "H")  # Increase count if roll == 1
  beta <- beta + (flips[i] != "H")    # Increase count otherwise
  posterior_means[i] <- alpha / (alpha + beta)  # Compute posterior mean
}

# Plot posterior mean convergence
plot(1:n_flips, posterior_means, type = "l", col = "blue", lwd = 2,
     xlab = "Number of Rolls", ylab = "Posterior Mean of p(rolling a 1)",
     main = "Bayesian Law of Large Numbers")
abline(h = 1/2, col = "red", lty = 2, lwd = 2)  # True probability reference line

library(ggplot2)

# Parameters
n <- 100
p <- 0.6

# Quantiles
q_low  <- qbinom(0.025, size = n, prob = p)
q_high <- qbinom(0.975, size = n, prob = p)

# PMF data frame
x_vals <- 0:n
df <- data.frame(
  x    = x_vals,
  prob = dbinom(x_vals, size = n, prob = p),
  region = case_when(
    x_vals < q_low               ~ "Lower Tail (2.5%)",
    x_vals > q_high              ~ "Upper Tail (2.5%)",
    TRUE                         ~ "Middle 95%"
  )
)

# Plot
ggplot(df, aes(x = x, y = prob, fill = region, color = region)) +
  geom_segment(aes(xend = x, yend = 0), linewidth = 0.8) +
  geom_point(size = 2) +
  geom_vline(xintercept = c(q_low, q_high), linetype = "dashed", color = "grey30") +
  annotate("label", x = q_low,  y = max(dbinom(x_vals, n, p)) * 0.9,
           label = paste0("q = ", q_low, "\n(p = 0.025)"),  size = 3.5, color = "#3b82f6") +
  annotate("label", x = q_high, y = max(dbinom(x_vals, n, p)) * 0.9,
           label = paste0("q = ", q_high, "\n(p = 0.975)"), size = 3.5, color = "#ef4444") +
  scale_color_manual(values = c(
    "Lower Tail (2.5%)" = "#3b82f6",
    "Middle 95%"        = "#22c55e",
    "Upper Tail (2.5%)" = "#ef4444"
  )) +
  scale_fill_manual(values = c(
    "Lower Tail (2.5%)" = "#3b82f6",
    "Middle 95%"        = "#22c55e",
    "Upper Tail (2.5%)" = "#ef4444"
  )) +
  labs(
    title    = "Binomial PMF with 95% Central Region",
    subtitle = paste0("X ~ Binomial(n = 100, p = 0.6)  |  95% interval: [", q_low, ", ", q_high, "]"),
    x        = "x",
    y        = "P(X = x)",
    fill     = NULL,
    color    = NULL
  ) +
  xlim(30, 90) +           # zoom in to relevant range
  theme_minimal(base_size = 14) +
  theme(legend.position = "top")

library(ggplot2)
library(dplyr)

# ============================================================
# SETUP: Human Height ~ Normal(mean = 170cm, sd = 10cm)
# ============================================================
mu <- 170   # mean height in cm
sigma <- 10 # standard deviation

# ============================================================
# 1. rnorm() — SIMULATE data (like sampling real people)
# ============================================================
# "If I randomly sampled 1000 people, what heights would I see?"

set.seed(42)
sample_heights <- rnorm(n = 1000, mean = mu, sd = sigma)

df_sample <- data.frame(height = sample_heights)

p1 <- ggplot(df_sample, aes(x = height)) +
  geom_histogram(aes(y = after_stat(density)), bins = 40,
                 fill = "#93c5fd", color = "white") +
  geom_density(color = "#1d4ed8", linewidth = 1) +
  stat_function(fun = dnorm, args = list(mean = mu, sd = sigma),
                color = "#ef4444", linewidth = 1, linetype = "dashed") +
  labs(title = "rnorm() — Simulated Sample (n = 1,000)",
       subtitle = "Blue = observed density | Red dashed = true Normal(170, 10)",
       x = "Height (cm)", y = "Density") +
  theme_minimal(base_size = 13)

# ============================================================
# 2. dnorm() — DENSITY (height of the curve at a point)
# ============================================================
# "How relatively likely is it to be exactly 170cm vs 190cm?"
# Note: this is NOT a probability — it's a density value.
# For continuous distributions, P(X = x) = 0 exactly.

x_vals <- seq(130, 210, length.out = 500)
df_density <- data.frame(
  x       = x_vals,
  density = dnorm(x_vals, mean = mu, sd = sigma)
)

# Points to highlight
highlight <- data.frame(
  x     = c(170, 180, 190),
  y     = dnorm(c(170, 180, 190), mean = mu, sd = sigma),
  label = c("dnorm(170) = 0.0399",
            "dnorm(180) = 0.0242",
            "dnorm(190) = 0.0054")
)

p2 <- ggplot(df_density, aes(x = x, y = density)) +
  geom_line(color = "#1d4ed8", linewidth = 1.2) +
  geom_point(data = highlight, aes(x = x, y = y), color = "#ef4444", size = 3) +
  geom_label(data = highlight, aes(x = x, y = y, label = label),
             nudge_y = 0.003, size = 3.2) +
  labs(title = "dnorm() — Probability Density Function (PDF)",
       subtitle = "Height of the curve — higher = more 'likely' region, but not a probability",
       x = "Height (cm)", y = "Density") +
  theme_minimal(base_size = 13)

# ============================================================
# 3. pnorm() — CDF & P-VALUES
# ============================================================
# "What is the probability of being shorter than X cm?"
# This is where p-values live — they are areas under the curve.

# --- One-tailed: P(Height < 185) ---
cutoff <- 185
p_lower <- pnorm(cutoff, mean = mu, sd = sigma, lower.tail = TRUE)
# "What % of people are shorter than 185cm?"

df_shade_lower <- df_density %>% filter(x <= cutoff)

p3 <- ggplot(df_density, aes(x = x, y = density)) +
  geom_area(data = df_shade_lower, aes(x = x, y = density),
            fill = "#3b82f6", alpha = 0.4) +
  geom_line(color = "#1d4ed8", linewidth = 1.2) +
  geom_vline(xintercept = cutoff, linetype = "dashed", color = "#ef4444") +
  annotate("label", x = 158, y = 0.025,
           label = paste0("P(X < 185) = ", round(p_lower, 4),
                          "\npnorm(185, 170, 10, lower.tail=T)"),
           size = 3.5, color = "#1d4ed8") +
  labs(title = "pnorm() — CDF: P(Height < 185cm)",
       subtitle = "Shaded area = probability = 0.9332 → ~93% of people are shorter than 185cm",
       x = "Height (cm)", y = "Density") +
  theme_minimal(base_size = 13)

# ============================================================
# 4. qnorm() — QUANTILES (inverse of pnorm)
# ============================================================
# "What height separates the top 5% from everyone else?"
# This is the CRITICAL VALUE concept in hypothesis testing.

q_95   <- qnorm(0.95,  mean = mu, sd = sigma)  # one-tailed 5%
q_025  <- qnorm(0.025, mean = mu, sd = sigma)  # two-tailed 2.5% lower
q_975  <- qnorm(0.975, mean = mu, sd = sigma)  # two-tailed 2.5% upper

df_shade_95   <- df_density %>% filter(x >= q_95)
df_shade_low  <- df_density %>% filter(x <= q_025)
df_shade_high <- df_density %>% filter(x >= q_975)

p5 <- ggplot(df_density, aes(x = x, y = density)) +
  geom_area(data = df_shade_low,  aes(x = x, y = density), fill = "#f97316", alpha = 0.5) +
  geom_area(data = df_shade_high, aes(x = x, y = density), fill = "#f97316", alpha = 0.5) +
  geom_line(color = "#1d4ed8", linewidth = 1.2) +
  geom_vline(xintercept = c(q_025, q_975), linetype = "dashed", color = "#f97316") +
  annotate("label", x = q_025, y = 0.015,
           label = paste0("qnorm(0.025)\n= ", round(q_025, 1), "cm"),
           size = 3.2, color = "#7c2d12") +
  annotate("label", x = q_975, y = 0.015,
           label = paste0("qnorm(0.975)\n= ", round(q_975, 1), "cm"),
           size = 3.2, color = "#7c2d12") +
  annotate("text", x = 170, y = 0.022,
           label = "Middle 95%\n(±1.96 SDs)", size = 3.5, color = "#1d4ed8") +
  labs(title = "qnorm() — Quantiles / Critical Values",
       subtitle = "qnorm(0.025) & qnorm(0.975) give the 95% interval — the origin of the ±1.96 rule!",
       x = "Height (cm)", y = "Density") +
  theme_minimal(base_size = 13)